MODULE SURFWS_INIT_ML_MOD
CONTAINS

SUBROUTINE SURFWS_INIT_ML(KIDIA, KFDIA, KLON, KLEVSN, NCL, PMU0,PSDOR,              &  ! Input
                     &  PTSOIL, PTSKIN, &
                     &  ZDSNTOT, ZSNDEPTH,                  &       
                     &  ZSNPERT,                                           &  ! Input
                     &  ZDSNR,PTSN, PRSN, PSSN, PWSN,PASN,                 &  ! Input
                     &  PTSNWS,PSSNWS,PRSNWS,PWSNWS,                       & ! Output
                     &  PTSNTOP, PTSNBOTTOM, PTSNMIDDLE,                   & ! Output
                     &  PRSNMAX, ZSADEPTH, KLEVSNA, KLEVMID, ZACTDEPTH,    & ! Output
                     &  PTMINCL,PRMINCL,                                   & ! Output
                     &  PTCONSTAVG, PTCONSTSTD,                            & ! Output
                     &  PRCONSTAVG, PRCONSTSTD, PRSNTOP,                   & ! Output
                     &  YDCST, YDSOIL )

USE PARKIND1 , ONLY : JPIM, JPRB
USE YOMHOOK  , ONLY : LHOOK, DR_HOOK, JPHOOK
USE YOS_CST  , ONLY : TCST
USE YOS_SOIL , ONLY : TSOIL

USE ABORT_SURF_MOD

! (C) Copyright 2017- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.

!**** *SURFWS_INIT_ML* - Snow warm start multi-layer 
!     PURPOSE.
!     --------
!          THIS ROUTINE CONTROLS THE WARM START OF MULTI-LAYER SCHEME IN
!          FC MODE

!**   INTERFACE.
!     ----------
!          *SURFWS_CTL* IS CALLED FROM *SRFSN_DRIVER*.

!     PARAMETER   DESCRIPTION                                    UNITS
!     ---------   -----------                                    -----

!     INPUT PARAMETERS (INTEGER):
!    *KIDIA*      START POINT


!     INPUT PARAMETERS (REAL):
!    *PTMST*      TIME STEP                                      S

!     INPUT PARAMETERS (LOGICAL):
!    *LDLAND*     LAND/SEA MASK (TRUE/FALSE)

!     INPUT PARAMETERS AT T-1 OR CONSTANT IN TIME (REAL):
!    *PWLM1M*     SKIN RESERVOIR WATER CONTENT                   kg/m**2


!     OUTPUT PARAMETERS AT T+1 (UNFILTERED,REAL):
!    *PWL*        SKIN RESERVOIR WATER CONTENT                   kg/m**2


!     OUTPUT PARAMETERS (DIAGNOSTIC):

!    *PDHIIS*     Diagnostic array for interception layer (see module yomcdh)

!     METHOD.
!     -------
!          

!     EXTERNALS.
!     ----------
!          NONE.

!     REFERENCE.
!     ----------
!          

!     Modifications:
!     Original   G. Arduini      ECMWF     28/07/2017

!     ------------------------------------------------------------------

IMPLICIT NONE

! Declaration of arguments 
INTEGER(KIND=JPIM), INTENT(IN) :: KIDIA
INTEGER(KIND=JPIM), INTENT(IN) :: KFDIA
INTEGER(KIND=JPIM), INTENT(IN) :: KLON
INTEGER(KIND=JPIM), INTENT(IN) :: KLEVSN
INTEGER(KIND=JPIM), INTENT(IN) :: NCL

REAL(KIND=JPRB), INTENT(IN)    :: ZDSNTOT(:), PTSN(:), PRSN(:), PSSN(:), PWSN(:)
REAL(KIND=JPRB), INTENT(IN)    :: ZSNDEPTH(:,:)
REAL(KIND=JPRB), INTENT(IN)    :: ZDSNR(:,:)
REAL(KIND=JPRB), INTENT(IN)    :: PSDOR(:)

REAL(KIND=JPRB), INTENT(IN)    :: PMU0(:)
REAL(KIND=JPRB), INTENT(IN)    :: PTSOIL(:)
REAL(KIND=JPRB), INTENT(IN)    :: PTSKIN(:)
REAL(KIND=JPRB), INTENT(IN)    :: PASN(:)

TYPE(TCST)     , INTENT(IN) :: YDCST
TYPE(TSOIL)    , INTENT(IN) :: YDSOIL

! Output Variables
REAL(KIND=JPRB), INTENT(OUT)    :: PTSNWS(:,:)
REAL(KIND=JPRB), INTENT(OUT)    :: PSSNWS(:,:)
REAL(KIND=JPRB), INTENT(OUT)    :: PRSNWS(:,:)
REAL(KIND=JPRB), INTENT(OUT)    :: PWSNWS(:,:)

REAL(KIND=JPRB), INTENT(OUT)    :: PTSNTOP(:), PTSNBOTTOM(:), PTSNMIDDLE(:)
REAL(KIND=JPRB), INTENT(OUT)    :: PRSNTOP(:)

REAL(KIND=JPRB), INTENT(OUT)    :: PRSNMAX(:)
REAL(KIND=JPRB), INTENT(OUT)    :: ZSADEPTH(:)

! Active number of layers:
INTEGER(KIND=JPIM),INTENT(INOUT)  :: KLEVSNA(:)

! Active depth layer (glaciers):
REAL(KIND=JPRB)    :: ZACTDEPTH(:)

! Cluster values:
INTEGER(KIND=JPIM),            INTENT(OUT)  :: KLEVMID(:)
INTEGER(KIND=JPIM),            INTENT(OUT)  :: PTMINCL(:), PRMINCL(:) 
REAL(KIND=JPRB),               INTENT(OUT)  :: PRCONSTAVG(:,:), PRCONSTSTD(:,:)
REAL(KIND=JPRB),               INTENT(OUT)  :: PTCONSTAVG(:,:), PTCONSTSTD(:,:)

REAL(KIND=JPRB)                 :: ZSNPERT

! Local variables:
INTEGER(KIND=JPIM)              :: JL,JK,KNACC
INTEGER(KIND=JPIM)              :: KCL
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZTDIST, ZRDIST
REAL(KIND=JPRB)                 :: ZFEATA, ZFEATB, ZFEATC 

! Look-up tables for warm-start
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZTCENTRADAY2, ZTCENTRBDAY2 ,ZTCENTRCDAY2 ,ZTCENTRADAY3 ,ZTCENTRBDAY3 ,ZTCENTRCDAY3 ,&
                                    ZTCENTRADAY4 ,ZTCENTRBDAY4 ,ZTCENTRCDAY4 ,ZTCENTRADAY5 ,ZTCENTRBDAY5 ,ZTCENTRCDAY5 ,&
                                    ZTCENTRADAY5M ,ZTCENTRBDAY5M ,ZTCENTRCDAY5M
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZTCENTRANIGHT2 ,ZTCENTRBNIGHT2 ,ZTCENTRCNIGHT2 ,ZTCENTRANIGHT3 ,&
                                    ZTCENTRBNIGHT3 ,ZTCENTRCNIGHT3 ,ZTCENTRANIGHT4 ,ZTCENTRBNIGHT4 ,&
                                    ZTCENTRCNIGHT4 ,ZTCENTRANIGHT5 ,ZTCENTRBNIGHT5 ,ZTCENTRCNIGHT5 ,&
                                    ZTCENTRANIGHT5M ,ZTCENTRBNIGHT5M ,ZTCENTRCNIGHT5M 
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZTCONSTAVGDAY2 ,ZTCONSTSTDDAY2 ,ZTCONSTAVGDAY3 ,ZTCONSTSTDDAY3 ,&
                                    ZTCONSTAVGDAY4 ,ZTCONSTSTDDAY4 ,ZTCONSTAVGDAY5 ,ZTCONSTSTDDAY5 ,&
                                    ZTCONSTAVGDAY5M ,ZTCONSTSTDDAY5M
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZTCONSTAVGNIGHT2 ,ZTCONSTSTDNIGHT2 ,ZTCONSTAVGNIGHT3 ,ZTCONSTSTDNIGHT3,&
                                    ZTCONSTAVGNIGHT4 , ZTCONSTSTDNIGHT4 ,ZTCONSTAVGNIGHT5 ,ZTCONSTSTDNIGHT5,&
                                    ZTCONSTAVGNIGHT5M ,ZTCONSTSTDNIGHT5M

REAL(KIND=JPRB), DIMENSION(NCL)  :: ZTCONSTAVG, ZRCONSTAVG, ZTCONSTSTD, ZRCONSTSTD
REAL(KIND=JPRB), DIMENSION(NCL)  ::  ZTCENTRA, ZTCENTRB, ZTCENTRC
REAL(KIND=JPRB), DIMENSION(NCL)  ::  ZRCENTRA, ZRCENTRB, ZRCENTRC


REAL(KIND=JPRB), DIMENSION(NCL)  :: ZTMLRIDAY2 ,ZTMLRADAY2 ,ZTMLRBDAY2 ,ZTMLRCDAY2 ,ZTMLRDDAY2 ,ZTMLREDAY2 
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZTMLRINIGHT2 ,ZTMLRANIGHT2 ,ZTMLRBNIGHT2 ,ZTMLRCNIGHT2 ,ZTMLRDNIGHT2 ,ZTMLRENIGHT2
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZTMLRIDAY3 ,ZTMLRADAY3 ,ZTMLRBDAY3 ,ZTMLRCDAY3 ,ZTMLRDDAY3 ,ZTMLREDAY3
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZTMLRINIGHT3 ,ZTMLRANIGHT3 ,ZTMLRBNIGHT3 ,ZTMLRCNIGHT3 ,ZTMLRDNIGHT3 ,ZTMLRENIGHT3
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZRCENTRA2 ,ZRCENTRB2 ,ZRCENTRC2 
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZRCENTRA3 ,ZRCENTRB3 ,ZRCENTRC3 
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZRCENTRA4 ,ZRCENTRB4 ,ZRCENTRC4
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZRCENTRA5 ,ZRCENTRB5 ,ZRCENTRC5
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZRCENTRA5M ,ZRCENTRB5M ,ZRCENTRC5M
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZRCONSTAVG2 ,ZRCONSTSTD2 ,ZRCONSTAVG3 ,ZRCONSTSTD3 ,&
                                    ZRCONSTAVG4 ,ZRCONSTSTD4 ,ZRCONSTAVG5 ,ZRCONSTSTD5 ,ZRCONSTAVG5M ,ZRCONSTSTD5M 
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZRMLRI2 ,ZRMLRA2 ,ZRMLRB2 ,ZRMLRC2 ,ZRMLRD2 ,ZRMLRE2 
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZRMLRI3 ,ZRMLRA3 ,ZRMLRB3 ,ZRMLRC3 ,ZRMLRD3 ,ZRMLRE3 
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZRMLRI4 ,ZRMLRA4 ,ZRMLRB4 ,ZRMLRC4 ,ZRMLRD4 ,ZRMLRE4 
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZRMLRI5 ,ZRMLRA5 ,ZRMLRB5 ,ZRMLRC5 ,ZRMLRD5 ,ZRMLRE5 
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZRMLRI5M ,ZRMLRA5M ,ZRMLRB5M ,ZRMLRC5M ,ZRMLRD5M ,ZRMLRE5M

REAL(KIND=JPRB), DIMENSION(NCL)  ::  ZTMLRA, ZTMLRB, ZTMLRC, ZTMLRD, ZTMLRE, ZTMLRI
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZRMLRI, ZRMLRA, ZRMLRB, ZRMLRC, ZRMLRD, ZRMLRE
REAL(KIND=JPRB)                  :: ZP1, ZP2, ZP3, ZP4, ZP5
REAL(KIND=JPRB)                  :: ZSDORTHR
REAL(KIND=JPRB)                  :: ZEPSILON

REAL(KIND=JPHOOK)                  :: ZHOOK_HANDLE


!    -----------------------------------------------------------------
IF (LHOOK) CALL DR_HOOK('SURFWS_INIT_ML_MOD:SURFWS_INIT_ML',0,ZHOOK_HANDLE)

!    -----------------------------------------------------------------

ASSOCIATE(RTT=>YDCST%RTT,RLMLT=>YDCST%RLMLT, RPI=>YDCST%RPI,    &
        & RLWCSWEA=>YDSOIL%RLWCSWEA, RLWCSWEB=>YDSOIL%RLWCSWEB, &
        & RLWCSWEC=>YDSOIL%RLWCSWEC, RTEMPAMP=>YDSOIL%RTEMPAMP, &
        & RDSNMAX=>YDSOIL%RDSNMAX, RHOMINSND=>YDSOIL%RHOMINSND, &
        & RHOMAXSN_NEW=>YDSOIL%RHOMAXSN_NEW, RDAT=>YDSOIL%RDAT  )

! mid-point first soil level
ZSADEPTH(:)=0.5_JPRB*RDAT(1)
ZSDORTHR=50._JPRB

ZEPSILON  = 10E4*EPSILON(ZEPSILON)

! 0.1 Define centroids 
! A: soT-skt B: soT-snT ; C: rsn/rsnmin
    ZTCENTRADAY2(1:NCL)=(/ 10.665,-0.746,-11.857,7.005,20.408,-4.571,-0.414,1.902,12.724,7.379,24.801,-9.194,-15.782,2.555,5.719,-5.452,3.056,17.634,-7.019,-2.571,-2.65,3.952,15.21,0.762,8.964,9.363,4.337  /)
    ZTCENTRBDAY2(1:NCL)=(/ 4.781,-0.267,-1.156,2.776,8.314,-1.646,0.063,2.046,5.265,3.294,9.803,-1.012,-1.084,0.515,2.159,-0.613,1.134,7.235,-1.633,-0.1,-0.92,1.695,6.044,0.423,9.546,3.372,3.996  /)
    ZTCENTRCDAY2(1:NCL)=(/ 2.494,2.595,4.17,5.094,2.082,2.899,5.544,2.275,2.339,2.231,1.962,4.18,4.105,2.342,2.275,4.885,4.896,2.155,3.299,5.529,2.843,2.294,2.205,2.509,1.935,2.248,2.103  /)
    ZTCENTRADAY3(1:NCL)=(/ 3.681,18.971,-7.308,11.437,0.007,16.422,24.398,6.981,-0.621,-13.271,13.776,-4.043,5.222,-9.71,33.778,9.633,2.111,-1.882,1.674,28.233,-6.001,2.332,7.409,-2.863,16.404,21.543,9.165  /)
    ZTCENTRBDAY3(1:NCL)=(/ 1.192,7.303,-1.917,4.373,-0.11,6.485,9.431,2.607,-0.017,-1.34,5.264,-1.859,1.918,-1.341,12.739,2.014,2.218,-1.048,0.498,10.614,-1.13,0.585,2.138,-0.295,4.092,8.052,3.804  /)
    ZTCENTRCDAY3(1:NCL)=(/ 2.402,2.155,3.088,2.423,2.713,2.235,2.124,2.307,5.504,4.062,2.335,2.788,2.35,3.974,1.913,2.228,2.292,2.84,2.427,2.007,4.273,4.91,5.157,5.258,1.887,2.144,2.445  /)
    ZTCENTRADAY4(1:NCL)=(/ 10.56,0.149,19.497,-9.183,3.922,35.092,16.334,3.376,-1.685,25.12,12.907,1.618,8.435,-1.946,6.902,17.105,-3.95,9.611,14.728,12.546,5.064,28.846,19.603,-6.319,-0.144,-12.722,22.192  /)
    ZTCENTRBDAY4(1:NCL)=(/ 3.897,-0.062,7.189,-1.702,0.519,11.943,3.312,1.252,-0.882,8.776,2.498,0.892,3.351,-0.026,2.114,6.206,-1.564,1.795,5.447,4.808,1.871,9.923,4.285,-1.773,0.064,-1.04,7.703  /)
    ZTCENTRCDAY4(1:NCL)=(/ 2.596,2.81,2.3,3.588,5.09,2.175,2.107,2.485,2.878,2.26,2.292,2.541,2.568,5.43,2.593,2.339,3.21,2.714,2.489,2.522,2.576,2.273,1.983,3.4,5.285,4.341,2.261  /)
    ZTCENTRADAY5(1:NCL)=(/ 32.408,6.341,-4.106,18.849,11.891,-11.135,-0.284,0.342,28.252,24.528,9.695,-2.07,7.203,-7.373,39.383,-1.44,3.577,14.854,3.385,5.561,7.991,14.081,21.734,10.979,18.98,1.741,16.374  /)
    ZTCENTRBDAY5(1:NCL)=(/ 8.96,3.454,-1.153,6.721,4.02,-1.178,0.145,0.195,7.992,7.344,3.599,0.2,0.071,-1.379,10.499,-0.702,1.948,2.654,0.051,1.695,2.396,4.975,6.354,1.03,4.506,1.251,5.356  /)
    ZTCENTRCDAY5(1:NCL)=(/ 2.403,2.48,3.326,2.347,2.571,3.956,5.006,2.974,2.425,2.408,2.58,5.37,4.161,3.54,2.379,3.426,2.587,2.853,3.927,2.73,2.732,2.477,2.421,3.58,2.545,2.665,2.448  /)
    ZTCENTRADAY5M(1:NCL)=(/ 4.785,19.234,0.792,8.587,26.497,13.488,6.71,12.668,-0.841,37.772,16.331,22.685,2.224,0.08,31.056,3.944,12.227,24.207,3.064,1.24,15.649,19.373,9.904,-1.57,7.331,10.082,5.897  /)
    ZTCENTRBDAY5M(1:NCL)=(/ 0.591,1.374,0.773,2.902,5.902,5.384,4.648,0.698,-0.028,7.553,4.797,5.595,0.503,0.49,6.097,3.673,3.204,1.952,1.285,2.971,1.096,5.13,0.77,0.411,0.659,5.089,2.142  /)
    ZTCENTRCDAY5M(1:NCL)=(/ 4.257,3.848,3.403,3.139,3.028,2.929,3.068,4.158,3.774,2.891,3.089,3.018,4.565,4.746,3.064,3.206,3.231,3.753,3.189,3.483,3.932,3.062,4.097,5.224,4.193,2.922,3.071  /)


    ZTCENTRANIGHT2(1:NCL)=(/ -0.651,14.093,6.256,9.776,19.373,-8.883,3.688,26.036,0.561,7.497,-5.54,17.692,0.838,5.362,22.849,3.867,10.815,30.659,15.924,3.715,21.54,-2.544,12.537,8.269,-2.698,2.156,7.078  /)
    ZTCENTRBNIGHT2(1:NCL)=(/ -0.186,6.389,3.488,4.084,8.645,-1.0,1.561,10.975,0.084,2.799,-1.574,7.641,0.511,2.239,9.833,3.08,5.212,11.767,6.944,0.616,7.849,-0.782,5.435,4.158,-0.997,1.358,2.534  /)
    ZTCENTRCNIGHT2(1:NCL)=(/ 2.432,2.2,2.179,2.264,2.111,4.34,2.241,1.945,5.193,2.196,3.266,2.151,2.307,2.302,2.009,2.201,2.388,1.949,2.178,5.111,1.856,4.485,2.223,2.336,2.441,2.373,4.791  /)
    ZTCENTRANIGHT3(1:NCL)=(/ 6.557,21.885,-0.557,12.786,30.201,16.116,9.468,2.275,17.978,-3.906,24.217,7.965,0.38,10.992,3.268,-8.062,-4.655,14.491,35.26,19.872,3.779,0.939,5.143,9.434,26.98,-2.078,6.393  /)
    ZTCENTRBNIGHT3(1:NCL)=(/ 2.757,8.877,-0.337,5.145,11.772,6.608,2.722,1.255,7.294,-0.915,9.685,3.386,-0.055,4.611,0.156,-1.426,-1.748,5.851,13.022,8.135,1.688,0.524,2.325,3.996,10.567,-1.156,0.624  /)
    ZTCENTRCNIGHT3(1:NCL)=(/ 2.32,2.118,2.533,2.349,2.025,2.172,4.207,2.368,2.171,4.513,2.06,2.294,5.187,2.278,5.071,3.454,2.597,2.18,1.883,2.179,2.31,2.495,2.267,2.25,1.976,3.035,4.949  /)
    ZTCENTRANIGHT4(1:NCL)=(/ 10.346,18.492,1.675,30.468,-3.631,11.878,22.46,7.479,42.471,16.697,-7.42,0.183,27.332,7.652,0.331,13.44,3.599,34.273,6.228,-1.455,4.756,11.882,9.075,20.353,15.044,24.729,3.274  /)
    ZTCENTRBNIGHT4(1:NCL)=(/ 3.969,6.637,0.823,10.492,-1.499,4.469,7.888,0.521,14.675,6.084,-1.503,-0.123,9.37,2.939,0.001,4.912,0.012,11.747,2.403,-1.125,1.895,2.52,3.261,7.353,5.421,8.676,1.346  /)
    ZTCENTRCNIGHT4(1:NCL)=(/ 2.451,2.285,2.456,2.129,3.237,2.427,2.236,4.481,1.999,2.276,3.316,2.533,2.169,2.463,5.218,2.43,4.848,2.106,2.452,2.669,2.469,3.239,2.6,2.239,2.418,2.235,2.42  /)
    ZTCENTRANIGHT5(1:NCL)=(/ 3.062,20.714,12.436,30.858,1.084,21.748,17.985,6.939,27.569,14.264,5.76,11.602,-5.275,10.587,8.819,3.14,34.898,14.982,24.912,8.65,0.465,4.99,41.511,16.229,18.683,-1.105,23.734  /)
    ZTCENTRBNIGHT5(1:NCL)=(/ 1.318,5.38,4.062,8.381,0.49,7.329,4.37,2.356,7.392,4.859,0.114,1.183,-1.27,3.606,2.901,-0.189,9.466,2.604,8.142,0.324,0.134,1.924,10.709,5.419,6.479,-0.929,5.933  /)
    ZTCENTRCNIGHT5(1:NCL)=(/ 2.542,2.471,2.516,2.336,2.672,2.287,2.674,2.561,2.394,2.424,4.304,3.558,3.697,2.497,2.583,4.339,2.348,3.061,2.337,4.267,4.969,2.505,2.301,2.406,2.316,3.132,2.434  /)
    ZTCENTRANIGHT5M(1:NCL)=(/ 8.979,19.759,0.575,29.292,15.061,14.129,6.453,40.803,15.995,10.028,23.28,2.944,3.806,25.493,10.616,7.739,5.169,12.397,33.645,1.44,6.382,21.858,26.95,18.062,21.286,12.538,17.699  /)
    ZTCENTRBNIGHT5M(1:NCL)=(/ 2.895,4.031,0.39,6.091,0.847,5.357,1.434,8.001,3.665,0.703,4.381,0.51,2.021,6.223,4.952,0.553,0.531,0.712,6.638,1.028,3.973,1.329,2.213,1.076,6.15,3.106,5.587  /)
    ZTCENTRCNIGHT5M(1:NCL)=(/ 3.097,3.375,4.654,3.04,4.053,2.851,3.108,2.855,3.303,4.199,3.333,4.552,2.971,2.913,2.867,4.519,4.497,4.146,2.957,3.127,2.969,3.854,3.697,3.948,2.822,3.271,2.905  /)


    ZTCONSTAVGDAY2(1:NCL)=(/ 4.053,0.873,-0.582,0.971,4.289,0.211,0.089,1.126,3.992,2.834,4.422,-0.625,-0.163,1.565,2.561,-0.059,0.935,4.022,-0.83,0.348,0.56,1.943,4.319,0.968,2.163,3.363,1.598  /)
    ZTCONSTSTDDAY2(1:NCL)=(/ 5.221,4.79,4.168,3.353,5.198,7.54,1.786,4.04,5.082,6.548,5.31,5.612,0.238,5.195,5.124,3.532,3.517,5.148,8.292,3.082,6.992,5.517,5.358,4.348,4.493,5.21,10.517  /)
    ZTCONSTAVGDAY3(1:NCL)=(/ 1.514,3.288,-1.483,2.978,0.562,3.338,3.753,2.541,0.032,-3.655,3.24,0.38,1.997,-2.318,3.914,3.128,0.664,0.53,0.587,3.856,-0.131,0.632,1.903,0.326,3.661,3.115,2.615  /)
    ZTCONSTSTDDAY3(1:NCL)=(/ 3.592,3.929,8.596,3.794,2.701,3.888,3.997,3.904,1.381,12.54,3.945,4.247,4.521,10.087,4.003,5.043,2.484,4.114,2.591,4.035,4.814,2.773,3.701,2.44,4.285,3.916,3.609  /)
    ZTCONSTAVGDAY4(1:NCL)=(/ 3.226,2.965,3.582,2.05,1.107,3.981,8.556,3.07,1.257,3.729,7.768,1.918,2.915,-0.066,5.013,3.719,0.234,8.105,3.597,3.298,3.306,3.973,7.663,1.725,0.145,1.475,4.03  /)
    ZTCONSTSTDDAY4(1:NCL)=(/ 4.591,7.93,3.857,4.382,3.25,3.245,7.857,6.983,5.359,3.551,8.707,5.796,4.258,0.173,7.249,4.054,2.71,8.987,4.029,3.887,6.317,3.179,7.319,18.775,2.04,3.857,3.882  /)
    ZTCONSTAVGDAY5(1:NCL)=(/ 4.481,0.507,-0.138,3.475,3.514,0.637,0.112,2.99,4.692,4.341,2.927,-0.052,6.375,0.072,4.95,1.855,1.76,7.815,4.916,4.722,4.467,3.387,4.416,7.049,6.319,1.844,3.821  /)
    ZTCONSTSTDDAY5(1:NCL)=(/ 2.493,5.065,1.455,2.944,4.206,1.753,1.926,7.789,3.894,2.881,4.658,0.781,7.116,1.766,1.906,6.553,6.253,5.788,7.905,6.731,5.814,3.517,3.073,6.346,4.195,6.194,3.689  /)
    ZTCONSTAVGDAY5M(1:NCL)=(/ 2.357,4.923,0.195,0.788,3.467,0.208,-1.324,4.882,0.533,4.068,1.609,2.87,0.928,-0.086,3.909,-1.59,1.709,4.662,0.728,-0.857,5.011,2.346,4.094,-0.092,3.482,-0.653,0.765  /)
    ZTCONSTSTDDAY5M(1:NCL)=(/ 4.144,2.206,3.196,3.147,1.685,2.034,3.418,2.899,3.605,1.211,2.338,1.929,3.311,0.703,1.412,4.494,2.94,1.649,4.368,6.933,2.7,2.591,3.401,0.428,3.902,1.992,3.92  /)
    ZTCONSTAVGNIGHT2(1:NCL)=(/ 0.18,3.361,2.589,3.044,3.859,-1.387,1.037,4.404,0.54,2.315,-0.084,3.632,0.369,1.741,3.776,1.773,3.292,4.645,3.604,1.243,4.265,-0.025,3.222,3.128,0.18,0.709,2.204  /)
    ZTCONSTSTDNIGHT2(1:NCL)=(/ 2.1,4.471,5.182,4.841,4.816,7.204,3.661,5.067,3.566,4.067,7.609,4.664,2.559,3.605,4.903,4.661,4.382,5.007,4.494,4.164,4.923,3.934,5.792,5.1,4.661,3.755,4.325  /)
    ZTCONSTAVGNIGHT3(1:NCL)=(/ 2.22,3.46,0.316,2.91,4.055,3.179,2.993,0.599,3.625,0.166,3.629,2.611,0.227,3.035,0.737,-2.982,-1.896,3.21,2.849,3.197,1.154,0.281,1.687,3.02,3.714,-0.005,2.114  /)
    ZTCONSTSTDNIGHT3(1:NCL)=(/ 3.546,3.971,2.306,3.746,4.077,3.825,4.132,2.285,3.905,2.321,4.006,4.036,2.0,3.732,3.051,11.839,9.679,3.828,4.036,3.869,2.939,1.878,3.269,3.74,4.064,4.294,4.201  /)
    ZTCONSTAVGNIGHT4(1:NCL)=(/ 3.201,3.663,1.107,4.285,0.455,3.149,3.869,5.2,4.462,3.528,1.42,1.528,4.028,3.081,0.173,3.667,3.074,3.836,2.548,1.35,2.443,4.682,3.177,3.609,3.645,3.926,1.626  /)
    ZTCONSTSTDNIGHT4(1:NCL)=(/ 3.548,3.63,4.099,3.513,2.602,3.452,3.726,6.385,3.645,3.669,3.523,5.479,3.534,3.965,1.88,3.488,6.468,3.446,3.959,5.01,4.437,4.894,3.939,3.571,3.585,3.58,4.487  /)
    ZTCONSTAVGNIGHT5(1:NCL)=(/ 1.704,5.005,3.86,4.804,1.622,3.983,5.51,3.042,4.92,3.597,6.646,6.787,0.408,3.486,3.409,4.552,4.967,6.232,4.114,6.731,0.357,2.043,5.654,3.914,3.923,1.053,5.074  /)
    ZTCONSTSTDNIGHT5(1:NCL)=(/ 4.418,2.964,2.863,2.562,5.618,3.183,3.336,3.905,2.5,2.728,7.205,4.653,1.942,3.01,3.494,8.171,2.44,3.719,2.745,6.063,3.078,3.904,2.489,3.528,3.302,6.092,2.514  /)
    ZTCONSTAVGNIGHT5M(1:NCL)=(/ 2.238,4.264,0.069,4.216,6.736,1.569,3.272,4.495,3.941,7.124,4.431,2.693,0.624,3.6,0.741,7.253,5.684,7.046,4.34,0.191,-0.264,5.769,5.23,6.391,2.904,3.515,2.504  /)
    ZTCONSTSTDNIGHT5M(1:NCL)=(/ 3.049,1.977,1.688,1.594,2.457,2.398,4.353,1.087,2.662,4.156,1.384,4.935,3.684,1.425,2.522,5.193,5.876,3.253,1.154,3.623,2.497,1.559,1.162,2.094,1.519,2.697,1.958  /)




! 0.3 Define constants for multi-linear regression of temp 1st layer (only ! dsn<0.20):
ZTMLRIDAY2(1:NCL)=(/2273.776,231.424,-1140.277,1077.753,-2605.152,4792.979,3877.053,393.788,&
                  &-9534.744,-11849.607,748.751,213.009,621.934,-296.89,2569.689,-870.507,&
                  &642.106,775.133,18059.695,218.558,-31549.314,-2667.837,1528.997,&
                  &-1178.686,176.003,2903.15,-2404.29 /)
ZTMLRADAY2(1:NCL)=(/-21.912,-22.981,-18.611,-3.775,19.878,-23.781,-28.105,-6.132,-39.31,&
                  &9.863,-42.609,-27.523,-7.325,-4.131,-21.878,-15.574,7.981,-6.876,5.878,&
                  &-0.976,7.58,-53.716,-14.317,-3.881,0.0,-19.631,-22.49 /)
ZTMLRBDAY2(1:NCL)=(/0.003,-0.013,0.002,0.003,0.001,0.0,-0.008,-0.001,0.003,0.001,-0.003,&
                  &0.003,0.011,-0.002,0.002,0.005,0.008,0.0,0.001,-0.0,0.001,-0.004,0.001,&
                  &0.002,-0.0,-0.004,0.001 /)
ZTMLRCDAY2(1:NCL)=(/2.63,21.656,23.678,-3.807,-0.0,-15.106,0.118,3.81,105.791,79.674,37.359,&
                  &21.693,1.89,7.147,0.852,18.914,-10.844,1.569,-136.832,0.277,225.808,73.195,&
                  &0.27,12.696,0.0,-0.986,41.659 /)
ZTMLRDDAY2(1:NCL)=(/0.046,0.043,0.043,0.007,-0.035,0.052,0.052,0.011,0.081,-0.018,0.079,0.061,&
                  &0.016,0.007,0.045,0.038,-0.019,0.013,-0.014,0.003,-0.015,0.101,0.034,&
                  &0.008,-0.0,0.036,0.041 /)
ZTMLREDAY2(1:NCL)=(/-0.002,-0.037,-0.041,0.01,0.001,0.031,0.003,-0.004,-0.192,-0.148,-0.066,&
                  &-0.037,-0.0,-0.01,0.002,-0.033,0.023,-0.0,0.254,0.0,-0.413,-0.133,0.002,&
                  &-0.021,0.001,0.004,-0.076/)

ZTMLRINIGHT2(1:NCL)=(/25.363,-4915.228,-272.045,-10385.244,-749.825,934.962,-1360.875,246.995,&
                    &-1131.312,1656.874,1945.898,-4666.717,3152.294,1125.078,15549.324,&
                    &-8619.04,-27203.797,769.829,-540.542,521.518,6036.785,127.023,&
                    &6869.104,-1328.559,885.202,-368.013,7971.363 /)
ZTMLRANIGHT2(1:NCL)=(/-21.853,-3.088,-43.653,-29.03,-16.206,1.009,4.797,6.269,-13.664,&
                    &-2.686,24.144,-11.438,-3.761,9.435,-10.072,5.866,-22.161,1.538,&
                    &-27.429,5.321,-30.082,-8.239,-7.267,-32.86,2.993,-41.127,-7.185 /)
ZTMLRBNIGHT2(1:NCL)=(/0.0,0.008,-0.007,0.002,-0.0,0.003,0.009,0.003,0.001,0.004,-0.019,&
                    &0.012,-0.001,0.002,-0.001,0.001,0.0,0.002,-0.002,-0.0,0.006,0.002,&
                    &-0.0,-0.008,-0.001,-0.004,0.0 /)
ZTMLRCNIGHT2(1:NCL)=(/21.915,40.216,45.865,103.003,22.934,-7.261,7.449,-6.425,20.528,&
                    &-9.486,-36.946,44.819,-19.352,-15.056,-103.585,59.556,220.627,-6.459,&
                    &31.715,-7.285,-18.251,7.928,-44.393,43.387,-8.759,43.923,-51.684 /)
ZTMLRDNIGHT2(1:NCL)=(/0.042,0.006,0.082,0.06,0.03,-0.003,-0.013,-0.014,0.029,0.007,-0.05,&
                    &0.026,0.007,-0.023,0.019,-0.01,0.045,-0.004,0.052,-0.013,0.064,0.015,&
                    &0.017,0.061,-0.007,0.078,0.014 /)
ZTMLRENIGHT2(1:NCL)=(/-0.039,-0.073,-0.083,-0.187,-0.041,0.017,-0.012,0.014,-0.035,0.02,0.074,&
                    &-0.081,0.039,0.03,0.192,-0.11,-0.403,0.016,-0.057,0.016,0.037,-0.012,&
                    &0.085,-0.078,0.02,-0.079,0.098/)

ZTMLRIDAY3(1:NCL)=(/-873.735,16851.732,1769.369,-1208.98,595.977,892.861,9584.742,&
                  &515.852,1914.724,-6049.516,-1895.017,-389.405,1705.265,1272.757,&
                  &29055.08,947.23,1665.529,5029.054,-333.738,24368.355,-1959.943,&
                  &1115.958,55.659,-13399.817,857.495,-1397.23,5746.104 /)
ZTMLRADAY3(1:NCL)=(/10.619,-15.266,-10.879,32.929,-25.726,14.589,-17.776,13.11,-27.152,&
                  &-42.891,14.106,-13.758,33.093,-5.929,19.702,40.24,-0.844,-20.965,&
                  &14.592,18.751,-18.933,-42.283,-39.847,7.549,21.17,-4.047,-9.189 /)
ZTMLRBDAY3(1:NCL)=(/0.0,0.0,0.0,0.001,-0.005,-0.008,-0.006,0.002,-0.004,-0.004,0.0,-0.001,&
                  &-0.003,0.004,0.002,-0.008,0.001,0.003,-0.01,-0.001,0.001,-0.008,&
                  &-0.003,0.001,-0.013,-0.001,0.0 /)
ZTMLRCDAY3(1:NCL)=(/-3.02,-109.68,-2.124,-22.655,21.564,-19.85,-53.218,-15.717,12.086,&
                  &84.146,-0.053,17.464,-44.619,-3.168,-232.189,-45.785,-11.335,-17.445,&
                  &-10.411,-196.785,32.479,33.39,39.859,90.671,-25.283,16.025,-33.291 /)
ZTMLRDDAY3(1:NCL)=(/-0.019,0.029,0.021,-0.059,0.049,-0.03,0.034,-0.027,0.052,0.087,-0.025,&
                  &0.026,-0.06,0.012,-0.039,-0.077,0.002,0.041,-0.033,-0.036,0.038,0.081,&
                  &0.075,-0.014,-0.046,0.008,0.018 /)
ZTMLREDAY3(1:NCL)=(/0.007,0.206,0.006,0.041,-0.038,0.041,0.101,0.033,-0.019,-0.153,0.002,&
                  &-0.031,0.083,0.008,0.431,0.088,0.024,0.036,0.025,0.365,-0.058,-0.06,&
                  &-0.072,-0.162,0.052,-0.029,0.064 /)

ZTMLRINIGHT3(1:27)=(/-233.194,1311.565,554.687,-170.099,685.766,5874.629,7978.278,&
                    &-8229.753,626.97,-163.629,1671.225,1383.694,22954.838,-6401.198,&
                    &1016.932,15437.012,465.208,5006.659,231532.438,-4.443,1691.159,768.856,&
                    &-1382.435,-7084.328,7818.812,1345.252,-3449.839 /)
ZTMLRANIGHT3(1:NCL)=(/-22.918,17.157,-22.351,-3.505,15.358,-58.272,-53.205,-3.579,31.338,&
                     &-5.707,15.969,23.453,28.833,-54.761,23.377,-32.093,-3.417,17.97,&
                     &11.031,-13.668,10.341,7.027,11.562,4.886,-1.373,19.383,32.732 /)
ZTMLRBNIGHT3(1:NCL)=(/-0.008,-0.002,-0.006,-0.002,-0.007,0.003,-0.008,-0.003,-0.002,&
                     &-0.007,0.009,-0.004,-0.004,-0.003,-0.006,-0.003,-0.007,-0.001,&
                     &0.017,-0.006,0.022,0.007,0.0,-0.006,-0.009,0.0,0.007 /)
ZTMLRCNIGHT3(1:NCL)=(/25.269,-25.646,16.266,5.886,-19.287,11.638,-8.116,66.062,-32.666,&
                    &7.912,-27.288,-32.353,-194.61,100.024,-29.553,-83.072,0.738,-53.752,&
                    &-1712.245,14.486,-22.205,-11.763,0.009,49.625,-55.516,-28.058,-2.126 /)
ZTMLRDNIGHT3(1:NCL)=(/0.043,-0.034,0.047,0.007,-0.029,0.115,0.104,0.006,-0.064,0.011,-0.032,&
                    &-0.045,-0.059,0.106,-0.045,0.061,0.007,-0.032,-0.021,0.026,-0.021,&
                    &-0.014,-0.021,-0.011,0.002,-0.038,-0.072 /)
ZTMLRENIGHT3(1:NCL)=(/-0.045,0.051,-0.028,-0.01,0.038,-0.019,0.018,-0.121,0.063,-0.013,&
                    &0.055,0.063,0.361,-0.182,0.058,0.158,0.0,0.099,3.15,-0.025,0.045,0.024,&
                    &0.001,-0.09,0.105,0.055,0.008 /)

! 0.4 Define centroids for Density: no distinction day and night...
! A: soT-skt B: soT-snT ; C: rsn/rsnmin
    ZRCENTRA2(1:NCL)=(/ -1334.4003,39.1089,-278.9727,52.0789,-278.5716,-225.1497,-160.8523,-111.9612,-87.4052,-101.3092,-383.7298,-1452.0217,-14.646,&
            & -347.6173,-880.4038,-1900.1066,-317.7257,-121.4379,-182.8199,-145.8385,135.3359,-17.3195,-71.4526,-18.1462,-1986.8328,-1717.4885,-200.5988 /)
    ZRCENTRB2(1:NCL)=(/ -93.1618,-511.9772,112.7408,-307.3933,-36.7917,303.3502,-347.5862,-319.4463,-73.7488,-266.8603,-320.346,-61.1226,-98.4689,&
            & -75.558,-35.2212,-81.9165,-42.9894,-116.3698,75.5675,-229.2397,-75.1291,-121.2313,-436.9606,-58.2537,-100.0915,-83.4104,71.8058 /)
    ZRCENTRC2(1:NCL)=(/ -114.608,356.9427,-89.1381,244.3195,-46.3497,-198.4637,264.4161,233.614,83.9798,186.5994,235.6509,-76.3997,93.8912,&
            & 40.7618,-44.5638,-101.6433,29.9802,79.395,-48.5811,207.4799,-99.3011,101.7144,278.4697,94.8146,-119.1965,-102.4522,-58.6891 /)
    ZRCENTRA3(1:NCL)=(/ 25.4653,42.5357,-350.1476,-131.4225,-907.8261,-70.0909,-9.3308,-45.5561,-173.2549,-301.3501,-1178.9126,-1349.7,-143.8232,&
            & -297.9534,-660.5016,-361.4285,-751.5026,-311.9067,153.5239,-317.777,-1736.8997,-220.5645,-366.6283,-174.2635,-12.5807,-23.9191,-1835.1179 /)
    ZRCENTRB3(1:NCL)=(/ -735.9476,-589.9268,153.9806,-84.8514,-58.5689,-426.7906,-473.0437,-264.4092,-450.7872,-32.5541,1614.6577,-220.9118,-406.3354,&
            & 68.2563,-66.1744,-23.1239,-141.1931,-197.6385,-96.3175,56.5325,-203.8742,-81.022,-67.2143,-500.8691,-400.9148,-549.0667,-153.5632 /)
    ZRCENTRC3(1:NCL)=(/ 461.6404,428.3682,-56.5207,127.2309,-76.8245,291.3167,336.7247,117.6744,336.9456,104.5965,-979.4467,-269.9912,323.4098,&
            & -123.1248,-85.5404,-31.9709,-178.1698,107.325,-124.7519,33.3334,-250.4757,138.7651,85.6363,286.7885,303.1394,392.2979,-193.3137 /)
    ZRCENTRA4(1:NCL)=(/ 122.0858,102.0765,-23.4785,-901.1704,197.393,-534.9584,217.3448,-116.3348,-69.1806,-208.2838,28.6775,-197.3348,-27.4655,&
            & -39.1663,159.2888,-158.6045,-402.7171,-876.2857,-181.1423,4.9736,203.6913,-388.4806,162.4262,59.7774,-439.1377,53.5147,134.453 /)
    ZRCENTRB4(1:NCL)=(/ -1130.1204,-214.9015,-707.3914,-146.1729,-1085.5802,-260.8198,-1053.7144,-592.6353,-512.8043,-92.479,-619.1938,-594.7343,-617.0192,&
            & -42.8982,-894.8987,6.6923,-411.6457,1694.0453,-539.5295,-709.4647,-41.1601,-53.8569,-45.8718,-719.5661,-157.0266,-1052.7983,-963.4905 /)
    ZRCENTRC4(1:NCL)=(/ 865.7322,-277.9244,477.937,-195.528,860.1669,-318.3297,736.1193,386.66,442.0512,-21.0093,525.9864,341.3652,478.4144,&
            & -57.0301,674.5895,9.9413,217.6331,-1075.1743,353.5464,516.4683,-54.2287,-74.8558,-66.8295,504.3237,-204.7632,749.9641,768.9747 /)
    ZRCENTRA5(1:NCL)=(/ -78.2416,79.8943,478.4736,-947.4033,143.7977,770.4563,595.1889,118.1891,668.2419,222.212,175.8237,-446.0686,317.1188,&
            & 150.5984,-552.5686,367.8889,58.5015,54.8829,127.1084,344.5151,802.0078,-464.8672,334.8764,335.4649,431.8332,226.1047,-543.7169 /)
    ZRCENTRB5(1:NCL)=(/ -61.3632,-1115.7137,-2418.6194,-213.5683,-1477.851,-2600.4905,-2251.5745,-1509.5658,-67.043,-1273.4856,-1549.4794,-147.2789,-2053.3398,&
            & -1318.4346,-109.0775,-803.3618,-2018.9907,-1249.0404,-1231.2401,-1775.7405,-70.4858,-133.9943,-1899.3363,-1719.322,-2135.946,-1716.0746,-131.2933 /)
    ZRCENTRC5(1:NCL)=(/ -93.2451,665.2209,1579.4581,-264.6112,958.7708,1692.5409,1460.9047,975.3162,-92.9966,785.9039,1017.8736,-200.0265,1323.3951,&
            & 759.5018,-150.3069,418.6025,1304.0142,731.1124,770.5744,1165.0714,-94.3249,-181.9821,1263.5157,1095.4689,1411.6543,1126.5073,-176.8941 /)
    ZRCENTRA5M(1:NCL)=(/ 164.2728,478.4838,486.5152,-650.8392,510.9012,-480.5661,495.9621,629.5404,-375.5007,-106.2239,1141.8042,800.9686,-328.1226,&
            & 1030.1239,621.0803,421.8679,-543.0834,1653.4512,1122.3219,1147.2332,458.8235,389.8102,217.3004,6.7735,413.6297,429.7506,-390.074 /)
    ZRCENTRB5M(1:NCL)=(/ -89.0774,-1130.7178,-2015.2688,-974.7406,-1291.9874,-34.1895,-1810.7769,-1927.7893,-46.9212,-1621.4994,-2139.2485,-1465.7592,-49.2912,&
            & -1729.3983,-1645.7028,-1794.6561,-58.6204,-207.3452,-1897.2556,-1710.0441,-1868.3457,-154.4845,-111.0024,-1140.3229,-1916.4017,-1220.0568,-39.856 /)
    ZRCENTRC5M(1:NCL)=(/ -138.0297,536.3357,1166.8176,213.2206,665.4451,-54.2601,1028.2416,1083.8828,-74.1052,918.4942,1096.8785,584.6857,-76.9579,&
            & 782.8547,899.1383,1022.8773,-92.3444,-270.5975,923.0989,758.0795,1082.1656,-229.2333,-168.0047,553.5491,1127.9065,624.525,-62.1632 /)


    ZRCONSTAVG2(1:NCL)=(/ 0.131,0.127,0.119,0.127,0.129,0.126,0.127,0.125,0.127,0.13,0.124,0.123,0.117,&
            & 0.124,0.122,0.127,0.124,0.13,0.126,0.113,0.135,0.13,0.128,0.12,0.112,0.115,0.123 /)
    ZRCONSTSTD2(1:NCL)=(/ 0.112,0.103,0.094,0.109,0.116,0.115,0.112,0.112,0.102,0.112,0.097,0.104,0.096,&
            & 0.101,0.104,0.106,0.105,0.108,0.11,0.092,0.128,0.109,0.113,0.099,0.078,0.092,0.109 /)
    ZRCONSTAVG3(1:NCL)=(/ 3.718,2.551,2.276,2.311,3.141,3.364,2.53,3.256,3.051,1.493,2.406,1.857,3.24,&
            & 3.525,2.394,2.766,1.541,4.365,2.016,1.879,1.542,2.151,2.052,3.628,2.184,3.578,2.753 /)
    ZRCONSTSTD3(1:NCL)=(/ 11.869,11.256,4.831,9.649,6.88,12.126,7.7,7.36,12.404,3.873,4.886,17.179,14.876,&
            & 9.279,5.325,8.692,4.992,15.254,4.696,4.818,8.476,4.775,4.201,11.849,7.616,14.706,14.517 /)
    ZRCONSTAVG4(1:NCL)=(/ -0.314,9.861,0.238,3.575,-1.912,5.002,-0.754,-0.024,-1.872,7.366,-0.804,2.951,-0.751,&
            & 19.437,-1.187,-2.925,0.215,-2.482,-0.276,0.318,17.349,1.799,0.363,1.536,11.107,-0.732,-0.814 /)
    ZRCONSTSTD4(1:NCL)=(/ 19.762,42.535,22.045,29.644,8.325,39.854,16.734,20.57,5.358,45.166,19.788,31.513,17.641,&
            & 64.103,14.409,3.438,21.782,4.099,21.394,22.85,63.111,27.252,21.977,29.846,42.105,17.061,17.582 /)
    ZRCONSTAVG5(1:NCL)=(/ 1.359,0.547,0.919,5.549,0.988,0.412,0.851,1.28,3.19,0.895,1.016,5.857,1.196,&
            & 0.904,4.979,1.688,2.622,0.452,0.941,1.656,3.423,5.96,1.007,1.129,1.528,1.141,6.923 /)
    ZRCONSTSTD5(1:NCL)=(/ 1.715,4.341,2.314,28.69,2.853,2.242,3.136,1.898,7.858,6.807,4.264,24.862,2.196,&
            & 4.524,20.687,5.351,13.385,4.048,4.164,1.727,8.714,23.222,3.186,1.59,2.016,3.563,29.225 /)
    ZRCONSTAVG5M(1:NCL)=(/ 0.166,-1.562,-1.251,3.171,-1.108,-0.909,-1.306,-1.139,-0.847,-1.27,1.924,2.782,-1.4,&
            & 2.809,-1.217,-1.193,-1.226,0.146,2.895,3.468,-1.416,1.454,0.16,-1.637,-1.189,-1.509,-0.992 /)
    ZRCONSTSTD5M(1:NCL)=(/ 2.291,2.43,2.362,18.431,5.422,2.996,1.983,3.755,2.181,1.992,11.886,15.901,2.377,&
            & 16.243,3.309,2.575,1.414,20.209,15.961,18.174,2.535,7.708,3.522,2.562,2.004,4.159,3.408 /)


    ZRMLRI2(1:NCL)=(/ -1334.4003,39.1089,-278.9727,52.0789,-278.5716,-225.1497,-160.8523,-111.9612,-87.4052,-101.3092,-383.7298,-1452.0217,-14.646,&
            & -347.6173,-880.4038,-1900.1066,-317.7257,-121.4379,-182.8199,-145.8385,135.3359,-17.3195,-71.4526,-18.1462,-1986.8328,-1717.4885,-200.5988 /)
    ZRMLRA2(1:NCL)=(/ -93.1618,-511.9772,112.7408,-307.3933,-36.7917,303.3502,-347.5862,-319.4463,-73.7488,-266.8603,-320.346,-61.1226,-98.4689,&
            & -75.558,-35.2212,-81.9165,-42.9894,-116.3698,75.5675,-229.2397,-75.1291,-121.2313,-436.9606,-58.2537,-100.0915,-83.4104,71.8058 /)
    ZRMLRB2(1:NCL)=(/ -0.0695,0.9455,0.6927,1.4983,0.9503,0.7838,0.6082,0.5228,3.7694,1.2504,0.3362,-0.2449,1.4198,&
            & 0.3453,0.4675,0.8537,0.2048,0.8852,0.5074,4.6993,0.7006,1.4562,0.7441,2.0489,1.8832,0.6999,0.7845 /)
    ZRMLRC2(1:NCL)=(/ 5.9333,0.5787,1.0556,-0.0479,1.1481,0.4999,1.1792,1.0127,-0.6281,0.6736,2.0922,6.3621,-0.0436,&
            & 1.6911,3.632,7.4894,1.5347,0.6758,0.777,-0.5273,-0.0003,0.027,1.0529,-0.4116,6.8645,6.8346,0.7842 /)
    ZRMLRD2(1:NCL)=(/ 0.0013,-0.0003,0.0003,-0.0014,-0.0001,0.0003,0.0004,0.0006,-0.0078,-0.0011,0.001,0.0016,-0.0012,&
            & 0.001,0.0006,0.0001,0.0015,-0.0002,0.0007,-0.0104,0.0002,-0.0015,0.0,-0.0026,-0.001,0.0003,0.0 /)
    ZRMLRE2(1:NCL)=(/ -114.608,356.9427,-89.1381,244.3195,-46.3497,-198.4637,264.4161,233.614,83.9798,186.5994,235.6509,-76.3997,93.8912,&
            & 40.7618,-44.5638,-101.6433,29.9802,79.395,-48.5811,207.4799,-99.3011,101.7144,278.4697,94.8146,-119.1965,-102.4522,-58.6891 /)
    ZRMLRI3(1:NCL)=(/ 25.4653,42.5357,-350.1476,-131.4225,-907.8261,-70.0909,-9.3308,-45.5561,-173.2549,-301.3501,-1178.9126,-1349.7,-143.8232,&
            & -297.9534,-660.5016,-361.4285,-751.5026,-311.9067,153.5239,-317.777,-1736.8997,-220.5645,-366.6283,-174.2635,-12.5807,-23.9191,-1835.1179 /)
    ZRMLRA3(1:NCL)=(/ -735.9476,-589.9268,153.9806,-84.8514,-58.5689,-426.7906,-473.0437,-264.4092,-450.7872,-32.5541,1614.6577,-220.9118,-406.3354,&
            & 68.2563,-66.1744,-23.1239,-141.1931,-197.6385,-96.3175,56.5325,-203.8742,-81.022,-67.2143,-500.8691,-400.9148,-549.0667,-153.5632 /)
    ZRMLRB3(1:NCL)=(/ -0.3258,0.563,3.9868,2.4426,0.26,0.5979,0.6311,-0.5929,0.8009,3.0996,6.3318,-0.8356,1.3818,&
            & -0.5228,0.3938,0.4547,-1.6156,0.5506,0.6345,3.5991,-2.0758,2.9712,5.2841,-0.0161,0.746,0.9899,-0.6318 /)
    ZRMLRC3(1:NCL)=(/ 1.5445,0.7691,-0.1387,-0.116,3.9511,1.0236,0.8128,1.3331,1.2768,0.1949,0.0555,7.0111,0.8476,&
            & 1.816,2.9945,1.7093,5.0406,1.7063,-0.0002,-0.0161,9.2597,-0.0128,-0.1165,1.8934,0.6229,0.836,8.5533 /)
    ZRMLRD3(1:NCL)=(/ 0.002,0.0004,-0.0077,-0.0038,0.0008,0.0001,0.0001,0.0025,0.0001,-0.0053,-0.0141,0.0022,-0.0013,&
            & 0.0023,0.0006,0.0006,0.0033,0.0001,0.0002,-0.0068,0.0037,-0.0051,-0.012,0.0012,0.0001,-0.0008,0.0019 /)
    ZRMLRE3(1:NCL)=(/ 461.6404,428.3682,-56.5207,127.2309,-76.8245,291.3167,336.7247,117.6744,336.9456,104.5965,-979.4467,-269.9912,323.4098,&
            & -123.1248,-85.5404,-31.9709,-178.1698,107.325,-124.7519,33.3334,-250.4757,138.7651,85.6363,286.7885,303.1394,392.2979,-193.3137 /)
    ZRMLRI4(1:NCL)=(/ 122.0858,102.0765,-23.4785,-901.1704,197.393,-534.9584,217.3448,-116.3348,-69.1806,-208.2838,28.6775,-197.3348,-27.4655,&
            & -39.1663,159.2888,-158.6045,-402.7171,-876.2857,-181.1423,4.9736,203.6913,-388.4806,162.4262,59.7774,-439.1377,53.5147,134.453 /)
    ZRMLRA4(1:NCL)=(/ -1130.1204,-214.9015,-707.3914,-146.1729,-1085.5802,-260.8198,-1053.7144,-592.6353,-512.8043,-92.479,-619.1938,-594.7343,-617.0192,&
            & -42.8982,-894.8987,6.6923,-411.6457,1694.0453,-539.5295,-709.4647,-41.1601,-53.8569,-45.8718,-719.5661,-157.0266,-1052.7983,-963.4905 /)
    ZRMLRB4(1:NCL)=(/ 2.5723,-2.8199,0.4754,0.9539,2.5196,-1.5958,0.8704,0.342,2.1525,-0.5352,2.5465,-0.144,1.2648,&
            & -0.3989,1.2461,4.45,0.0552,3.1034,0.3544,1.3448,-0.434,0.1004,0.5158,0.3904,-0.8635,0.8139,2.5885 /)
    ZRMLRC4(1:NCL)=(/ 0.2128,3.0255,1.2539,3.9477,-0.2148,4.6587,0.5703,1.5526,0.2305,1.7284,-0.1753,2.1469,0.6698,&
            & 0.984,0.2903,-0.8257,2.5899,-0.0842,1.7077,0.7305,0.1301,2.1024,-0.3232,0.933,3.394,1.2102,-0.1367 /)
    ZRMLRD4(1:NCL)=(/ -0.0041,0.0048,0.0005,-0.0004,-0.0038,0.0032,-0.0001,0.0006,-0.0028,0.0022,-0.0039,0.0015,-0.001,&
            & 0.0019,-0.0007,-0.0097,0.001,-0.0055,0.0005,-0.0015,0.0019,0.0013,0.0008,0.0008,0.0024,-0.0001,-0.004 /)
    ZRMLRE4(1:NCL)=(/ 865.7322,-277.9244,477.937,-195.528,860.1669,-318.3297,736.1193,386.66,442.0512,-21.0093,525.9864,341.3652,478.4144,&
            & -57.0301,674.5895,9.9413,217.6331,-1075.1743,353.5464,516.4683,-54.2287,-74.8558,-66.8295,504.3237,-204.7632,749.9641,768.9747 /)
    ZRMLRI5(1:NCL)=(/ -78.2416,79.8943,478.4736,-947.4033,143.7977,770.4563,595.1889,118.1891,668.2419,222.212,175.8237,-446.0686,317.1188,&
            & 150.5984,-552.5686,367.8889,58.5015,54.8829,127.1084,344.5151,802.0078,-464.8672,334.8764,335.4649,431.8332,226.1047,-543.7169 /)
    ZRMLRA5(1:NCL)=(/ -61.3632,-1115.7137,-2418.6194,-213.5683,-1477.851,-2600.4905,-2251.5745,-1509.5658,-67.043,-1273.4856,-1549.4794,-147.2789,-2053.3398,&
            & -1318.4346,-109.0775,-803.3618,-2018.9907,-1249.0404,-1231.2401,-1775.7405,-70.4858,-133.9943,-1899.3363,-1719.322,-2135.946,-1716.0746,-131.2933 /)
    ZRMLRB5(1:NCL)=(/ 1.856,0.0652,2.4044,-0.0554,1.0751,2.2282,2.4132,1.7944,-1.3339,0.4019,1.2578,-1.8233,2.1692,&
            & -0.9206,-1.3642,-1.3662,0.0852,0.027,0.5736,1.7195,-1.7229,-1.5678,1.76,3.0594,1.5232,1.3761,-1.974 /)
    ZRMLRC5(1:NCL)=(/ 0.4096,1.8828,1.1244,4.7691,1.5511,0.3775,0.4626,1.4745,-1.0327,1.3805,1.4202,3.9259,1.3943,&
            & 2.3857,3.8504,1.0715,3.0812,2.243,1.5664,0.8513,-1.3221,3.7598,1.0699,0.4741,1.0941,1.4037,4.2574 /)
    ZRMLRD5(1:NCL)=(/ -0.0039,0.0005,-0.0044,0.0016,-0.0013,-0.0041,-0.0046,-0.003,0.0035,-0.0002,-0.0016,0.0039,-0.004,&
            & 0.0026,0.0034,0.0035,0.0011,0.0005,-0.0005,-0.0024,0.0042,0.0037,-0.0027,-0.0064,-0.0019,-0.0017,0.0042 /)
    ZRMLRE5(1:NCL)=(/ -93.2451,665.2209,1579.4581,-264.6112,958.7708,1692.5409,1460.9047,975.3162,-92.9966,785.9039,1017.8736,-200.0265,1323.3951,&
            & 759.5018,-150.3069,418.6025,1304.0142,731.1124,770.5744,1165.0714,-94.3249,-181.9821,1263.5157,1095.4689,1411.6543,1126.5073,-176.8941 /)
    ZRMLRI5M(1:NCL)=(/ 164.2728,478.4838,486.5152,-650.8392,510.9012,-480.5661,495.9621,629.5404,-375.5007,-106.2239,1141.8042,800.9686,-328.1226,&
            & 1030.1239,621.0803,421.8679,-543.0834,1653.4512,1122.3219,1147.2332,458.8235,389.8102,217.3004,6.7735,413.6297,429.7506,-390.074 /)
    ZRMLRA5M(1:NCL)=(/ -89.0774,-1130.7178,-2015.2688,-974.7406,-1291.9874,-34.1895,-1810.7769,-1927.7893,-46.9212,-1621.4994,-2139.2485,-1465.7592,-49.2912,&
            & -1729.3983,-1645.7028,-1794.6561,-58.6204,-207.3452,-1897.2556,-1710.0441,-1868.3457,-154.4845,-111.0024,-1140.3229,-1916.4017,-1220.0568,-39.856 /)
    ZRMLRB5M(1:NCL)=(/ -0.4118,-0.5071,1.3942,-2.9421,0.3626,3.3835,1.4226,0.571,2.2049,0.2784,-2.237,-3.0275,0.8641,&
            & -2.8014,0.4664,1.2754,6.0582,-1.2952,-2.7106,-2.7934,0.4291,-1.4963,-0.6344,-0.7916,1.8469,0.3031,0.9638 /)
    ZRMLRC5M(1:NCL)=(/ 0.829,1.0918,1.2287,6.8376,0.7197,0.9093,0.9186,1.0621,1.1798,3.3815,1.2714,2.1608,1.6524,&
            & 1.5024,0.7432,1.2444,0.1893,-3.9432,1.2858,1.0557,1.575,0.9551,0.8759,2.9339,1.0945,0.9477,1.7464 /)
    ZRMLRD5M(1:NCL)=(/ 0.0008,0.0011,-0.0023,0.0051,-0.0005,-0.0061,-0.0024,-0.0008,-0.0038,-0.0003,0.004,0.0053,-0.0014,&
            & 0.005,-0.0007,-0.0021,-0.0123,0.0033,0.0048,0.0049,-0.0006,0.0027,0.0012,0.0017,-0.0031,-0.0005,-0.0015 /)
    ZRMLRE5M(1:NCL)=(/ -138.0297,536.3357,1166.8176,213.2206,665.4451,-54.2601,1028.2416,1083.8828,-74.1052,918.4942,1096.8785,584.6857,-76.9579,&
            & 782.8547,899.1383,1022.8773,-92.3444,-270.5975,923.0989,758.0795,1082.1656,-229.2333,-168.0047,553.5491,1127.9065,624.525,-62.1632 /)



!*******************************************************************************
! 2. Start snow parametrizations
!    Here we start the snow temperature and density parametrisation:
!    both are simple exponential function. Temperature relaxes to first soil layer
!    at the bottom, while the density relaxes to mean density. 
!    The missing snow mass it is added to the bottom layer simply increasing the
!    density of this layer keeping the depth fixed.
!*******************************************************************************
  DO JL=KIDIA, KFDIA
    IF (PMU0(JL) > ZEPSILON ) THEN
      IF ( ZDSNTOT(JL) < 0.15_JPRB ) THEN 
        ZTCENTRA=ZTCENTRADAY2
        ZTCENTRB=ZTCENTRBDAY2
        ZTCENTRC=ZTCENTRCDAY2

        ZTCONSTAVG=ZTCONSTAVGDAY2
        ZTCONSTSTD=ZTCONSTSTDDAY2

        ZTMLRI=ZTMLRIDAY2
        ZTMLRA=ZTMLRADAY2
        ZTMLRB=ZTMLRBDAY2
        ZTMLRC=ZTMLRCDAY2
        ZTMLRD=ZTMLRDDAY2
        ZTMLRE=ZTMLREDAY2

      ELSEIF ( ZDSNTOT(JL) < 0.20_JPRB ) THEN 
        ZTCENTRA=ZTCENTRADAY3
        ZTCENTRB=ZTCENTRBDAY3
        ZTCENTRC=ZTCENTRCDAY3

        ZTCONSTAVG=ZTCONSTAVGDAY3
        ZTCONSTSTD=ZTCONSTSTDDAY3
        
        ZTMLRI=ZTMLRIDAY3
        ZTMLRA=ZTMLRADAY3
        ZTMLRB=ZTMLRBDAY3
        ZTMLRC=ZTMLRCDAY3
        ZTMLRD=ZTMLRDDAY3
        ZTMLRE=ZTMLREDAY3

      ELSEIF ( ZDSNTOT(JL) < 0.25_JPRB ) THEN 
        ZTCENTRA=ZTCENTRADAY4
        ZTCENTRB=ZTCENTRBDAY4
        ZTCENTRC=ZTCENTRCDAY4
        
        ZTCONSTAVG=ZTCONSTAVGDAY4
        ZTCONSTSTD=ZTCONSTSTDDAY4

      ELSEIF ( ZDSNTOT(JL) < 0.50_JPRB ) THEN 
        ZTCENTRA=ZTCENTRADAY5
        ZTCENTRB=ZTCENTRBDAY5
        ZTCENTRC=ZTCENTRCDAY5

        ZTCONSTAVG=ZTCONSTAVGDAY5
        ZTCONSTSTD=ZTCONSTSTDDAY5
        
      ELSEIF ( ZDSNTOT(JL) >= 0.50_JPRB ) THEN 
        ZTCENTRA=ZTCENTRADAY5M
        ZTCENTRB=ZTCENTRBDAY5M
        ZTCENTRC=ZTCENTRCDAY5M

        ZTCONSTAVG=ZTCONSTAVGDAY5M
        ZTCONSTSTD=ZTCONSTSTDDAY5M
      ENDIF
        
    ELSEIF (PMU0(JL)<=ZEPSILON) THEN
      IF ( ZDSNTOT(JL) < 0.15_JPRB ) THEN 
        ZTCENTRA=ZTCENTRANIGHT2
        ZTCENTRB=ZTCENTRBNIGHT2
        ZTCENTRC=ZTCENTRCNIGHT2

        ZTCONSTAVG=ZTCONSTAVGNIGHT2
        ZTCONSTSTD=ZTCONSTSTDNIGHT2

        ZTMLRI=ZTMLRINIGHT2
        ZTMLRA=ZTMLRANIGHT2
        ZTMLRB=ZTMLRBNIGHT2
        ZTMLRC=ZTMLRCNIGHT2
        ZTMLRD=ZTMLRDNIGHT2
        ZTMLRE=ZTMLRENIGHT2
        
      ELSEIF ( ZDSNTOT(JL) < 0.20_JPRB ) THEN 
        ZTCENTRA=ZTCENTRANIGHT3
        ZTCENTRB=ZTCENTRBNIGHT3
        ZTCENTRC=ZTCENTRCNIGHT3

        ZTCONSTAVG=ZTCONSTAVGNIGHT3
        ZTCONSTSTD=ZTCONSTSTDNIGHT3
        
        ZTMLRI=ZTMLRINIGHT3
        ZTMLRA=ZTMLRANIGHT3
        ZTMLRB=ZTMLRBNIGHT3
        ZTMLRC=ZTMLRCNIGHT3
        ZTMLRD=ZTMLRDNIGHT3
        ZTMLRE=ZTMLRENIGHT3

      ELSEIF ( ZDSNTOT(JL) < 0.25_JPRB ) THEN 
        ZTCENTRA=ZTCENTRANIGHT4
        ZTCENTRB=ZTCENTRBNIGHT4
        ZTCENTRC=ZTCENTRCNIGHT4

        ZTCONSTAVG=ZTCONSTAVGNIGHT4
        ZTCONSTSTD=ZTCONSTSTDNIGHT4
        
      ELSEIF ( ZDSNTOT(JL) < 0.50_JPRB ) THEN 
        ZTCENTRA=ZTCENTRANIGHT5
        ZTCENTRB=ZTCENTRBNIGHT5
        ZTCENTRC=ZTCENTRCNIGHT5

        ZTCONSTAVG=ZTCONSTAVGNIGHT5
        ZTCONSTSTD=ZTCONSTSTDNIGHT5
        
      ELSEIF ( ZDSNTOT(JL) >= 0.50_JPRB ) THEN 
        ZTCENTRA=ZTCENTRANIGHT5M
        ZTCENTRB=ZTCENTRBNIGHT5M
        ZTCENTRC=ZTCENTRCNIGHT5M

        ZTCONSTAVG=ZTCONSTAVGNIGHT5M
        ZTCONSTSTD=ZTCONSTSTDNIGHT5M
      ENDIF
    ENDIF

    PTCONSTAVG(JL,1:NCL)=ZTCONSTAVG(1:NCL)
    PTCONSTSTD(JL,1:NCL)=ZTCONSTSTD(1:NCL)

! Assign density depending on snow depth::
    IF ( ZDSNTOT(JL) < 0.15_JPRB ) THEN 
      ZRCENTRA=ZRCENTRA2
      ZRCENTRB=ZRCENTRB2
      ZRCENTRC=ZRCENTRC2

      ZRCONSTAVG=ZRCONSTAVG2
      ZRCONSTSTD=ZRCONSTSTD2

      ZRMLRI=ZRMLRI2
      ZRMLRA=ZRMLRA2
      ZRMLRB=ZRMLRB2
      ZRMLRC=ZRMLRC2
      ZRMLRD=ZRMLRD2
      ZRMLRE=ZRMLRE2
      
      KLEVMID(JL)=2_JPIM
    ELSEIF ( ZDSNTOT(JL) < 0.20_JPRB ) THEN 
      ZRCENTRA=ZRCENTRA3
      ZRCENTRB=ZRCENTRB3
      ZRCENTRC=ZRCENTRC3

      ZRCONSTAVG=ZRCONSTAVG3
      ZRCONSTSTD=ZRCONSTSTD3
      
      ZRMLRI=ZRMLRI3
      ZRMLRA=ZRMLRA3
      ZRMLRB=ZRMLRB3
      ZRMLRC=ZRMLRC3
      ZRMLRD=ZRMLRD3
      ZRMLRE=ZRMLRE3

      KLEVMID(JL)=2_JPIM

    ELSEIF ( ZDSNTOT(JL) < 0.25_JPRB ) THEN 
      ZRCENTRA=ZRCENTRA4
      ZRCENTRB=ZRCENTRB4
      ZRCENTRC=ZRCENTRC4

      ZRCONSTAVG=ZRCONSTAVG4
      ZRCONSTSTD=ZRCONSTSTD4
      
      ZRMLRI=ZRMLRI4
      ZRMLRA=ZRMLRA4
      ZRMLRB=ZRMLRB4
      ZRMLRC=ZRMLRC4
      ZRMLRD=ZRMLRD4
      ZRMLRE=ZRMLRE4
 
      KLEVMID(JL)=3_JPIM

    ELSEIF ( ZDSNTOT(JL) < 0.50_JPRB ) THEN 
      ZRCENTRA=ZRCENTRA5
      ZRCENTRB=ZRCENTRB5
      ZRCENTRC=ZRCENTRC5

      ZRCONSTAVG=ZRCONSTAVG5
      ZRCONSTSTD=ZRCONSTSTD5
      
      ZRMLRI=ZRMLRI5
      ZRMLRA=ZRMLRA5
      ZRMLRB=ZRMLRB5
      ZRMLRC=ZRMLRC5
      ZRMLRD=ZRMLRD5
      ZRMLRE=ZRMLRE5

      IF (PSDOR(JL)<ZSDORTHR)THEN
        KLEVMID(JL)=3_JPIM
      ELSEIF (KLEVSNA(JL)<KLEVSN .AND. KLEVSNA(JL)>2)THEN
        KLEVMID(JL)=KLEVSNA(JL)-1
      ELSEIF (KLEVSNA(JL)==2)THEN
        KLEVMID(JL)=KLEVSNA(JL)
      ELSE
        KLEVMID(JL)=3_JPIM
      ENDIF

    ELSEIF ( ZDSNTOT(JL) >= 0.50_JPRB ) THEN 
      ZRCENTRA=ZRCENTRA5M
      ZRCENTRB=ZRCENTRB5M
      ZRCENTRC=ZRCENTRC5M

      ZRCONSTAVG=ZRCONSTAVG5M
      ZRCONSTSTD=ZRCONSTSTD5M
      
      ZRMLRI=ZRMLRI5M
      ZRMLRA=ZRMLRA5M
      ZRMLRB=ZRMLRB5M
      ZRMLRC=ZRMLRC5M
      ZRMLRD=ZRMLRD5M
      ZRMLRE=ZRMLRE5M

      KLEVMID(JL)=MAX(KLEVSNA(JL)-1,1)
     
    ENDIF
! Final assignment:
!  RCENTRA(JL)=ZRCENTRA
!  RCENTRB(JL)=ZRCENTRB
!  RCENTRC(JL)=ZRCENTRC

   PRCONSTAVG(JL,1:NCL)=ZRCONSTAVG(1:NCL)
   PRCONSTSTD(JL,1:NCL)=ZRCONSTSTD(1:NCL)
   
!  RMLRI(JL)=ZRMLRI
!  RMLRA(JL)=ZRMLRA
!  RMLRB(JL)=ZRMLRB
!  RMLRC(JL)=ZRMLRC
!  RMLRD(JL)=ZRMLRD
!  RMLRE(JL)=ZRMLRE

!*****************************************
! 2.1 Find closest cluster centre using euclidean metric:
!     Features are common to temp and dens
   ZFEATA=PTSOIL(JL)-PTSKIN(JL)
   ZFEATB=PTSOIL(JL)-PTSN(JL)
   ZFEATC=PRSN(JL)/RHOMINSND

! 2.1.1 Temperature
   DO KCL=1,NCL
     ZTDIST(KCL)=SQRT( (ZFEATA-ZTCENTRA(KCL))**2_JPRB + (ZFEATB-ZTCENTRB(KCL))**2_JPRB + (ZFEATC-ZTCENTRC(KCL))**2_JPRB )
   ENDDO
   PTMINCL(JL)=MINLOC(ZTDIST,DIM=1)

! 2.1.2 Density:
   DO KCL=1,NCL
     ZRDIST(KCL)=SQRT( (ZFEATA-ZRCENTRA(KCL))**2_JPRB + (ZFEATB-ZRCENTRB(KCL))**2_JPRB + (ZFEATC-ZRCENTRC(KCL))**2_JPRB )
   ENDDO
   PRMINCL(JL)=MINLOC(ZRDIST,DIM=1)

!*****************************************
! 2.1 Initialize warm start (WS) variables
    IF ( PSSN(JL) < ZSNPERT ) THEN
        PTSNWS(JL, 1)          = PTSKIN(JL)
        PTSNWS(JL, 2:KLEVSN-1) = PTSN(JL)
        PTSNWS(JL, KLEVSN)     = PTSOIL(JL)

        PRSNWS(JL, 1:KLEVSN) = PRSN(JL)
        PRSNMAX(JL)          = RHOMAXSN_NEW

        PSSNWS(JL, 1)        = PSSN(JL)
        PWSNWS(JL, 1)        = PWSN(JL)
        PSSNWS(JL, 2:KLEVSN) = 0._JPRB
        PWSNWS(JL, 2:KLEVSN) = 0._JPRB

!----- Density
        ZP1=ZRMLRI(PRMINCL(JL))+ZRMLRA(PRMINCL(JL))*PASN(JL)
        ZP2=                    ZRMLRB(PRMINCL(JL))*PRSN(JL)
        ZP3=                    ZRMLRC(PRMINCL(JL))*PTSN(JL)
        ZP4=                    ZRMLRD(PRMINCL(JL))*(PRSN(JL))**2_JPRB
        ZP5=                    ZRMLRE(PRMINCL(JL))*(PASN(JL))**2_JPRB

        PRSNTOP(JL)=ZP1+ZP2+ZP3+ZP4+ZP5

        PTSNBOTTOM(JL)=PTSOIL(JL)
        PTSNMIDDLE(JL)=0.5*(PTSN(JL)+PTSOIL(JL))
        PTSNTOP(JL)   =PTSN(JL)
        IF ( ZDSNTOT(JL) < 0.20_JPRB ) THEN
          ZSADEPTH(JL)=0._JPRB
        ELSE
          ZSADEPTH(JL)=0.5_JPRB*RDAT(1)
          IF (PTSNTOP(JL)>=PTSN(JL)) THEN
            ZACTDEPTH(JL)=ZSNDEPTH(JL, KLEVSNA(JL)-1)
          ELSE
            ZACTDEPTH(JL)=ZDSNTOT(JL)
          ENDIF
        ENDIF

    ELSE ! Glacier ini
        KLEVSNA(JL)=KLEVSN !KSNACC+1

        PTSNWS(JL, 1)        = PTSKIN(JL)
        PTSNWS(JL, 2:KLEVSN) = PTSN(JL)

        PRSNMAX(JL)              = 300._JPRB
        PRSNWS(JL,1:KLEVSN)      = PRSNMAX(JL)
       !PRSNWS(JL,KSNACC:KLEVSN) = PRSNMAX(JL)
        PSSNWS(JL,1:KLEVSN)      = PRSNMAX(JL)*ZDSNR(JL,1:KLEVSN)
        PWSNWS(JL,1:KLEVSN)      = 0._JPRB
!----- Active depth
      PTSNTOP(JL)   = PTSKIN(JL)
      IF (PTSNTOP(JL)>=PTSN(JL)) THEN
        ZACTDEPTH(JL)=ZSNDEPTH(JL, KLEVSNA(JL))
      ELSE
        ZACTDEPTH(JL)=ZDSNTOT(JL)
      ENDIF
        PTSNMIDDLE(JL)=0.5_JPRB*(PTSN(JL)+PTSOIL(JL))
        PTSNBOTTOM(JL)=PTSOIL(JL)
!----- Density
      PRSNTOP(JL)=300._JPRB
    ENDIF
  ENDDO 
END ASSOCIATE

!    -----------------------------------------------------------------
IF (LHOOK) CALL DR_HOOK('SURFWS_INIT_ML_MOD:SURFWS_INIT_ML',1,ZHOOK_HANDLE)

END SUBROUTINE SURFWS_INIT_ML
END MODULE SURFWS_INIT_ML_MOD


